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We propose a new generalized-ensemble algorithm, which we refer to as the multibaric- 
multithermal Monte Carlo method. The multibaric-multithermal Monte Carlo simulations perform 
random walks widely both in volume space and in potential energy space. From only one simulation 
run, one can calculate isobaric-isothermal-ensemble averages at any pressure and any temperature. 
We test the effectiveness of this algorithm by applying it to the Lennard- Jones 12-6 potential system 
with 500 particles. It is found that a single simulation of the new method indeed gives accurate 
average quantities in isobaric-isothermal ensemble for a wide range of pressure and temperature. 
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Monte Carlo (MC) algorithm is one of the most widely used methods of computational physics. In order to realize 
desired statistical ensembles, corresponding MC techniques have been proposed [II 0> Oa S 13 • The first MC simulation 
was performed in the canonical ensemble in 1953 This method is called the Metropolis algorithm and widely used. 
The canonical probability distribution Pj^yt{E) for potential energy E is given by the product of the density of states 
n(E) and the Boltzmann weight factor e~ l3aE : 

Pnvt(E) = n{E)e-f 3oE , (1) 

where (3q is the inverse of the product of the Boltzmann constant ks and temperature To &t which simulations are 
performed. Since n(E) is a rapidly increasing function and the Boltzmann factor decreases exponentially, P^vt(E) 
is a bell-shaped distribution. 

The isobaric-isothermal MC simulation 2] is also extensively used. This is because most experiments are carried out 
under the constant pressure and constant temperature conditions. Both potential energy E and volume V fluctuate 
in this ensemble. The distribution P^pt(E, V) for E and V is given by 

P NPT (E,V)=n(E,V)e-' 3 ° H . (2) 

Here, the density of states n(E, V) is given as a function of both E and V, and H is the "enthalpy" : 

H = E + P Q V , (3) 

where Pq is the pressure at which simulations are performed. This ensemble has bell-shaped distributions in both E 
and V. 

Besides the above physical ensembles, it is now almost a default to simulate in artificial, generalized ensembles so 
that the multiple- minima problem, or the broken ergodicity problem, in complex systems can be overcome (for a recent 
review, see Ref. @). The multicanonical algorithm 0,0 is one of the most well known such methods in generalized 
ensemble. In multicanonical ensemble, a non-Boltzmann weight factor W mc (E) is used. This multicanonical weight 
factor is characterized by a flat probability distribution P mc (E): 



S ' Pmc(£) = n{E)W mc {E) = constant , (4) 

_ _ i 

and thus a free random walk is realized in the potential energy space. This enables the simulation to escape from any 
local-minimum-energy state and to sample the configurational space more widely than the conventional canonical MC 
algorithm. Another advantage is that one can obtain various canonical ensemble averages at any temperature from a 
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single simulation run by the reweighting techniques [9j- This method is now widely used in complex systems such as 
proteins and glasses However, it is difficult to compare the simulation conditions with experimental environments 
of constant pressure, since the simulations are performed in a fixed volume. 

In this Letter, we propose a new MC algorithm in which one can obtain various isobaric-isothermal ensembles 
from only one simulation. In other words, we introduce the idea of the multicanonical technique into the isobaric- 
isothermal ensemble MC method. We call this method the multibaric-multithermal algorithm. This MC simulation 
performs random walks in volume space as well as in potential energy space. As a result, this method has the 
following advantages: (1) It allows the simulation to escape from any local-minimum-energy state and to sample 
the configurational space more widely than the conventional isobaric-isothermal method. (2) One can obtain various 
isobaric-isothermal ensembles not only at any temperature, as in the multicanonical algorithm, but also at any 
pressure from only one simulation run. (3) One can control pressures and temperatures similarly to real experimental 
conditions. 

In the multibaric-multithermal ensemble, every state is sampled by a weight factor W m bt(E,V) = 
exp{— /3oH m ht(E, V)} (H m bt is referred to as the multibaric-multithermal enthalpy) so that a uniform distribution in 
both potential energy space and volume space is obtained: 

Pmbt(£, V) = n{E, V)W mht (E, V) = constant . (5) 

We call W m bt{E,V) the multibaric-multithermal weight factor. 

In order to perform the multibaric-multithermal MC simulation, we follow the conventional isobaric-isothermal MC 
techniques |2J. In this method, we perform Metropolis sampling on the scaled coordinates Si — L~ 1 r i (r^ are the real 
coordinates) and the volume V (here, the particles are placed in a cubic box of a side of size L = \/V). The trial 
moves of the scaled coordinates from Si to s'i and of the volume from V to V are generated by uniform random 
numbers. The enthalpy is accordingly changed from H(E(s^ N ' , V), V) to H'(E(s'( N \ V'), V) by these trial moves. 
The trial moves will be accepted with the probability 

acc(o -> n) = min(l, exp[-/? {# ' -H- Nk B T ln(V7V)}]) , (6) 

where N is the total number of particles in the system. 

Replacing H by -f/mbt, we can perform the multibaric-multithermal MC simulation. The trial moves of Si and V 
are generated in the same way as in the isobaric-isothermal MC simulation. The multibaric-multithermal enthalpy 
is changed from H m bt(E(s( N \ V), V) to H' mht {E{s'^ N \ V), V) by these trial moves. The trial moves will now be 
accepted with the probability 

acc(o - n) = min(l, exp[-p {H^ bt - H mht - Nk B T \n(V'/V)}]) . (7) 

The multibaric-multithermal probability distribution P m bt(^, V) is obtained by this scheme. 

In order to calculate the isobaric-isothermal ensemble average, we employ the reweighting techniques 0. The prob- 
ability distribution Pnpt(E, V;T, P) at any temperature T and any pressure P in the isobaric-isothermal ensemble 
is given by 

**M*. v; t. p) - T Y tE ' V,Wjt(E - V,e "" B+PV ' ■ («> 

I dV I dE P m u(E,V) W,-i,(E,V) e-'^ FV > 
The expectation value of a physical quantity A at T and P is estimated from 

< 



:A> NPT = J dV J dE A(E,V)P N pt(E,V;T,P) , 

^ AfT? T/\W-1 I J? Xl\„-0(E+PV) ^ . 



< A{E,V)W-l t {E,V)e^ E + pv ) > mbt 

<W-£ t (E,V)e-P( E + p v) > mbt ' U 

where < • • • > m bt is the multibaric-multithermal ensemble average. 

After having given the formalism of the multibaric-multithermal algorithm, let us now describe the process for 
determining the weight factor W m bt(E, V). This is obtained by the usual iteration of short simulations [lOL UlL Il2j|. 
The first simulation is carried out at To and Pq in the isobaric-isothermal ensemble. Namely, we use 

W£1(E, V) = exp{-p H^ t (E, V)} , (10) 
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where 



H™ t (E,V)=E + P V . (11) 



The i-th simulation is performed with the weig ht factor w£l t (E, V) and let P® U (E, V) be the obtained distribution. 
The (i + l)-th weight factor (E, V) is then given by 

W%£\E,V) = eM-PoH^iE^)} , (12) 



where 



Mi+Drw, ) H £l( E ,V) + k B T lnpW u (E,V) , for P» t (£, V) > , 

#mbt ( A > K ) = S /,n ( ld ) 

H$ ht (E,V), for P^(£?,V) = 0. 

For convenience, we make E and V discrete into bins and use histograms to calculate P® ht (E,V). We iterate this 
process until a reasonably flat distribution P^ ht {E,V) is obtained. After an optimal weight factor is determined, a 
long simulation is performed to sample the configurational space. 

We now present the results of our multibaric-multithcrmal simulation. We considered a Lennard- Jones 12-6 potential 
system. We used 500 particles (N = 500) in a cubic unit cell with periodic boundary conditions. The length and the 
energy are scaled in units of the Lennard- Jones diameter a and the minimum value of the potential e, respectively. 
We use an asterisk (*) for the reduced quantities such as the reduced length r* = r/a, the reduced temperature 
T* = kftT/e, the reduced pressure P* — Pa 3 /e, and the reduced number density p* = per 3 (p = N/V). 

We started the multibaric-multithermal weight factor determination of Eqs. (|12fl and (|13|l from a regular isobaric- 
isothcrmal simulation at Tg = 2.0 and P * = 3.0. These temperature and pressure are respectively higher than the 
critical temperature T! and the critical pressure P* [HQ 111 HE3 E3 . Recent reliable data are T* = 1.3207(4) 
and P* = 0.1288(5) [l^]. The cutoff radius r* was taken to be L*/2. A cut-off correction was added for the pressure 
and potential energy. In one MC sweep we made the trial moves of all particle coordinates and the volume (TV + 1 
trial moves altogether). For each trial move the Metropolis evaluation of Eq. 10 was made. In order to obtain a 
flat probability distribution P m ht(E,V), we carried out the MC simulations of 100,000 MC sweeps and iterated the 
process of Eqs. (|12fl and (|13fl . In the present case, it was required to make 12 iterations to get an optimal weight 
factor W m bt(E, V). We then performed a long multibaric-multithermal MC simulaton of 400,000 MC sweeps with 
this W mbt (E,V). 

For the purpose of comparisons of the new method with the conventional one, we also performed the conventional 
isobaric-isothermal MC simulations of 100,000 MC sweeps with 500 Lennard- Jones 12-6 potential particles at several 
sets of temperature and pressure. The temperature ranged from T* — 1.6 to 2.6 and the pressure from P* = 2.2 to 
3.8. 

In order to estimate the statistical accuracies, we performed these MC simulations from four different initial con- 
ditions in both multibaric-multithermal and isobaric-isothermal simulations. The error bars were calculated by the 
standard deviations from these different simulations. 

Figure 2] shows the probability distributions of E* /N and V*/N. Figure QJa) is the probability distribution 
Pnpt{E* /N, V*/N) from the isobaric-isothermal simulation first carried out in the process of Eqs. I|10fl and i|ll|) (i.e., 
T * = 2.0 and Pq = 3.0). It is a bell-shaped distribution. On the other hand, Fig. ^b) is the probability distribution 
Pmbt(E* /N, V*/N) from the multibaric-multithermal simulation finally performed. It shows a flat distribution, and 
the multibaric-multithermal MC simulation indeed sampled the configurational space in wider ranges of energy and 
volume than the conventional isobaric-isothermal MC simulation. 

Figure [5] shows the time series of E*/N. Figure Ufa) gives the results of the conventional isobaric-isothermal 
simulations at (T*,P*) = (1.6,3.0) and (2.4,3.0), while Figure Hfb) presents those of the multibaric-multithermal 
simulation. The potential energy fluctuates in narrow ranges in the conventional isobaric-isothermal MC simulations. 
They fluctuate only in the ranges of E* /N = —4.0 ~ —3.6 and E* /N — —5.1 ~ —4.7 at the higher temperature of 
T* = 2.4 and at the lower temperature of T* = 1.6, respectively. On the other hand, the multibaric-multithermal 
MC simulation performs a random walk that covers a wide energy range. 

A similar situation is observed in the time series of V* JN. In Fig. |3{ a) we show the results of the conventional 
isobaric-isothermal simulations at (T*,P*) — (2.0,2.2) and (2.0,3.8), while in Figure IUb) we give those of the 
multibaric-multithermal simulation. 

The volume fluctuations in the conventional isobaric-isothermal MC simulations are only in the range of V* /N — 
1.3 ~ 1.4 and V*/N = 1.5 - 1.6 at P* = 3.8 and at P* = 2.2, respectively. On the other hand, the multibaric- 
multithermal MC simulation performs a random walk that covers even a wider volume range. 
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We calculated the ensemble averages of potential energy per particle, < E* /N >npt, and density, < p* >npt, 
at various temperature and pressure values by the rcweighting techniques of Eq. @. They are shown in Fig. 01 
and in Fig. 03 respectively. The error bars are smaller than the plots for both cases. The agreement between the 
multibaric-multithermal data and isobaric-isothermal data are excellent in both < E* /N >npt and < p* >npt- 

The important point is that we can obtain any desired isobaric-isothermal distribution in wide temperature and 
pressure ranges (T* = 1.6 ~ 2.6, P* = 2.2 ~ 3.8) from a single simulation run by the multibaric-multithermal 
MC algorithm. This is an outstanding advantage over the conventional isobaric-isothermal MC algorithm, in which 
simulations have to be carried out separately at each temperature and pressure, because the reweighting techniques 
based on the isobaric-isothermal simulations can give correct results only for narrow ranges of temperature and 
pressure values. 

Figures 0] and [5] also show two equations of states of the Lennard- Jones 12-6 potential fluid. One is determined 
by Johnson et al. [Tg| and the other by Sun and Teja |20j |. These equations are determined by fitting procedure to 
the molecular simulation results. Our multibaric-multithermal simulation results agree very well with those of these 
equations. Investigating in more detail, however, the two equations give slightly different results. Most of our data 
lie in between them. 

In conclusion, we proposed a new MC algorithm that is based on multibaric-multithermal ensemble. We applied this 
method to the Lennard-Jones 12-6 potential system. The advantage of this method is that the simulation performs 
random walks in both potential energy space and volume space and sample the configurational space much more 
widely than the conventional isobaric-isothermal MC method. Therefore, one can obtain various isobaric-isothermal 
ensemble averages at any desired temperature and pressure from only one simulation run. This allows one to specify a 
pressure and to compare simulation conditions directly with those of real experiments. The multibaric-multithermal 
algorithm will thus be of great use for investigating a large variety of complex systems such as proteins, polymers, 
supercooled liquids, and glasses. It will be particularly useful for the study of, for example, pressure induced phase 
transitions. 

We would like to thank M. Mikami of National Institute of Advanced Industrial Science and Technology and U.H.E. 
Hansmann of Michigan Technological University for useful discussions at the early stage of the present work. 
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FIG. 1: (a) The probability distribution V NPT (E* /N, V/N) in the isobaric-isothermal simulation at (T* , P*) = (T * , P * ) = 
(2.0,3.0) and (b) the probability distribution P ulht (E* /N, V/N) in the multibaric-multithermal simulation. 
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FIG. 2: The time series of E* /N from (a) the conventional isobaric-isothermal MC simulations at (T*, P*) — (2.4, 3.0) and at 
(T*,P*) = (1.6,3.0) and (b) the multibaric-multithermal MC simulation. 
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FIG. 3: The time series of V* /N from (a) the conventional isobaric-isothermal MC simulations at (T*, P*) — (2.0, 2.2) and at 
(T*,P*) = (2.0,3.8) and (b) the multibaric-multithermal MC simulation. 
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FIG. 4: Average potential energy per particle < E* /N >npt at various temperature and pressure values. Filled circles: 
Multibaric-multithermal MC simulations. Open squares: Conventional isobaric-isothermal MC simulations. Solid line: Equa- 
tion of states calculated by Johnson et al. [19]. Broken line: Equation of states calculated by Sun and Teja [20]. 




FIG. 5: Average density < p* >npt at various temperature and pressure values. See the caption of Fig. 4 for details. 



